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^ Abstract 
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Among the efficient numerical methods based on atomistic models, the quasicontinuum (QC) 
method has attracted growing interest in recent years. The QC method was first developed for 
crystalline materials with Bravais lattice and was later extended to multilattices (Tadmor et 
al, 1999). Another existing numerical approach to modeling multilattices is homogenization. 
In the present paper we review the existing numerical methods for multilattices and propose 
, /^ i another concurrent macro-to-micro method in the homogenization framework. We give a unified 

(-H mathematical formulation of the new and the existing methods and show their equivalence. We 

then consider extensions of the proposed method to time-dependent problems and to random 
materials. 
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)0 1 Introduction 

■^ 

m 

In some applications of solid mechanics, such as modeling cracks, structural defects, or nanoelec- 

Q tromechanical systems, the classical continuum description is not suitable, and it is required to 

T-H utilize an atomistic description of materials. However, full atomistic simulations are prohibitively 

^7^ expensive, hence one needs to coarse-grain the problem. The quasicontinuum (QC) method [39] is 

K>" one of the most efficient methods of coarse-graining the atomistic statics. The idea behind QC is to 

K/ introduce a piecewise affine constraints for the atoms in regions with smooth deformation and use 

?-H the Cauchy-Born rule to define the energy of the corresponding groups of constrained atoms. To 

formulate the QC method for multilattice crystals one must account for relative shifts of Bravais 

lattices which the multilattice is comprised of [30] . 

The QC method is a multiscale method capable of coupling atomistic and continuum description 
of materials. It is intended to model an atomistic material in a continuum manner in the regions 
where the deformation is smooth and use the fully atomistic model only in the small neighborhood of 
defects, thus effectively reducing the degrees of freedom of the system. Originally, the QC method 
was developed for crystalline materials with a (single) Bravais lattice j39j and the convergence 
of a few variants of the method has been analyzed under some practical assumptions (see, e.g., 
[m [El [m EH Sni ED Ea ES]). The QC method is based on the so-called Cauchy-Born rule 
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(see, e.g., [9l [jQl [231 125]) which states that the energy of a certain volume of a material can be 
approximated through the deformation energy density, which is computed for a representative atom 
assuming that the neighboring atoms follow a uniform deformation. Later, QC was extended to 
multilattices [lO] (a multilattice is a union of a number of Bravais lattices) based on the improved 
Cauchy-Born rule |38] which accounts for relative shifts between the Bravais lattices. Examples of 
such materials include diamond cubic Si, HCP metals (stacking two simple hexagonal lattices with a 
shift vector) like Zr, ferroelectric materials, salts like Sodium Chloride, and intermet allies like NiAl. 
More recent developments of QC for multilattices also include adaptive choice of representative cell 
of multilattices [15j . It appears that no rigorous analysis is available so far for the multilattice QC. 

In the present work we propose a treatment of multilattices within the framework of numerical 
homogenization. Homogenization techniques for partial differential equations (PDEs) with mul- 
tiscale coefficients are known to be successful for obtaining effective equations with coefficients 
properly averaged out |8]. Finite element methods based on homogenization theory have been pio- 
neered by Babiiska [5] and have attracted growing attention these past few years (see [U [181 [HI [26] 
for textbooks or review papers). Following the ideas of |8j, we use homogenization techniques to 
describe the coarse-graining of multilattices and based on that, propose a macro-to-micro numeri- 
cal algorithm which we call the homogenized QC (HQC) method. We give a unified mathematical 
description and establish equivalence between the homogenized QC, the multilattice QC of |40j . 
and the finite element method applied to the continuously homogenized equations (see [H 123] and 
references therein for homogenization of atomistic media) . 

Despite the formal equivalence of these numerical methods for the multilattices, we find sev- 
eral benefits of the homogenization framework. First, in this framework the connection to the 
well-developed theory of continuum homogenization and related numerical methods becomes more 
apparent. This allows us to apply the numerical analysis techniques developed for continuum ho- 
mogenization [U [18] to the multilattice QC method (see the preprint [4J for an example of such 
application). Second, homogenization theory can be used to upscale the atomistic model in both, 
time and space, which makes it promising for modeling and especially analyzing zero temperature 
and finite temperature motion of atomistic materials |18l [^ [55] . In this work we demonstrate such 
an application of HQC to zero-temperature dynamics (Section p|. Also, homogenization can be 
applied to "stochastic" materials, atomistic counterparts of which include polymers [7] and glasses 
(we give an example of such application in Section uj). Last, for the finite temperature simulations, 
when materials are modeled with static atoms interacting with effective temperature-dependent 



potentials, the discrete homogenization (Section 4.3) may serve as a rigorous instrument to derive 



such potentials. We note that the idea of applying homogenization to atomistic media has appeared 
in the hterature before [3 [IH [HI [H [23] . 

The paper is organized as follows. We present the atomistic model in Section [2| and the 
quasicontinuum method in Section [3} In Section |4] we present the technique of homogenization 
applied to the atomistic equations. In Section [5] we present the HQC method — a concurrent macro- 
to-micro algorithm based on the discrete homogenization. The method is formulated in such a 
way that it allows for a straightforward extension to non-crystalline materials if the microstructure 
is known; an example of such extension is given in Section [7[ Section [6] is devoted to showing 
the equivalence of the following four methods: the HQC method, the multilattice QC, the finite 
element method applied to continuously homogenized equations, and the QC method applied to 
the discretely homogenized equations. In Section [8] we apply the proposed macro-to-micro method 
to a long-wave unsteady evolution of a ID multilattice crystal. Concluding remarks are given in 



Section [9j The commonly used notations are collected in the appendix. 

2 Problem Formulation 

The focus of the present study is on correct treatment of atomistic materials with spatially oscil- 
lating or inhomogeneous local properties. 

2.1 Equations of Equilibrium 

We describe the formulation of the problem of finding an equilibrium of an atomistic material in 
the periodic setting. We consider the periodic boundary conditions for simplicity, in order to avoid 
difficulties arising from presence of the boundary of the atomistic material. Nevertheless, it should 
be noted that the numerical method and the algorithm proposed in the present work can be applied 
to Dirichlet, Neumann, or other boundary conditions. 

2.1.1 Deformation 

Consider an atomistic material occupying a region f] = (0,1] in its reference (i.e., undeformed) 
configuration and extended periodically outside of O. The set of positions of atoms in the reference 
configurations are 

m— 1 



M = nn (j {eZ'^ + epa), 



a=0 

where pa S [0, l)*^ is a shift vector of a-th species of atoms in the reference configuration; in total 
we have m species of atoms. We assume that pa / P/3 for a / /? and, for convenience, po = 0. 

We collect these shift vectors into the set V := {pa '■ a = 0, . . . ,m — 1}. Thus, if we denote a 
Bravais lattice in fi as 

c = nn eZ"', 

then we can write A4 = C + eV. We will call M a multilattice. 

When the material experiences a deformation the atom positions become x + u{x), where u{x) 
is the displacement. We assume that u{x) is periodic, i.e., u{x + a) = u{x) for all a G Z . The space 
of all periodic displacements is denoted as Uperi-M)- Since we only consider the systems invariant 
w.r.t. translation in space, we will also need the space of displacements with zero average, U^{M.) 



(see Appendix A.l for the precise definitions) 



2.1.2 Interaction 

We assume the most general form of interaction (finite-range, multi-body) between atoms. For 
each atom x G A^ we introduce its "interaction neighborhood" — a set of vectors TZe{x) such 
that {x + er : r G TZe{x)} are the atoms that x interacts with. The energy of an atom x G Ai 



is denoted as Ve(L'7^^(a;)u(x); x), where D-ji^u = {Dru)reTie (see Appendix A. 3) is a collection of 
discrete directional derivatives of u corresponding to the set of neighbors TZe (these notations were 
first introduced in [28j ) . The discrete derivative in direction r of n evaluated at x E A^ is defined 
as Dru{x) := "i^ er)-u(x) ^ r^^^ needed properties and definitions of discrete directional derivatives 



can be found in Appendix ] A. 2 and more details on discrete directional derivatives — in Appendix 
Thus, the interaction energy of the displacement u is given by the interaction potential V^ as 

where {g)s denotes the average value of a function g defined on a discrete set S. 

The subscript e in V^ and TZ^ indicates that these objects depend non-smoothly on x: indeed, 
the interaction energy and the interaction neighborhood may depend on the species of atoms a 
for X £ C + epoi- For instance, we can consider a Lennard- Jones potentials with atom-dependent 
parameters: 

y(I)^.n;x)= 5:..,.,.(-2(g^)-V(g^)-^^), (2.1) 

relic 

where Sx,x+er and ix,x+er are, respectively, the strength and the equilibrium distance of interaction 
of atoms X and x + er. 

We assume that the interaction neighborhood TZ^ix+epa), and the interaction potential Ve{», x+ 
epa), for X G £ depends only on a, the particular species of atoms, but does not depend on x; we 
therefore write TZ^ix + epa) ='■ Tle^a and Ve{».,x + epa) ='■ V^^a- Then, we can use the following form 
of the energy: 

m— 1 



E{u) = (— X] ^^(^Tl,{x+ep„)u{x + epa); X + epa 
0=0 
^ m— 1 

Y,yUDTl.,Mx + epa))) , (2.2) 

fit I Xt- Lj 



where we used a more verbose notation for averaging of a function g defined on a discrete set S, 
{g)s ='■ {g{x))xes- This expression for the energy will be used to write down the energy of the 



multilattice QC method in a familiar way (see (3.5)). 



The latter form of the energy will be used whenever we need periodicity of interaction (for 
the Cauchy-Born rule and for the homogenization) . However we stress that the formulation of 
the proposed numerical method (Section ^ requires neither periodicity of Ai nor periodicity of 
interaction. 

2.1.3 External Force 

The potential energy of the external force / = /(x) is 

-F{u) = -{f,u)M, 

where by {w,v)m := {w ■ v)m we denote a scalar product oi w,v E Upcr{-M.)- (To be precise, it 
is an inner product on U^{M.) and a semi-inner product on ^pcr(A^).) The forces / = f{x) are 
applied as "dead loads", i.e., they are independent of actual atom positions x + u. For the problem 
to be well-posed, the sum of all forces per period is assumed to be zero, i.e., (/)x = 0. 



2.1.4 Equation of Equilibrium 

We denote the total potential energy of the atomistic system as 

n(n) = E{u)-F{u). 

A displacement u £ U^{M) is a stable equilibrium if it is a local minimizer of 11, which implies 
that It is a critical point of 11: 

{m{u),v)M-=^^^iu + tv)l^^ = yv£U#{M). (2.3) 

We assume that the function n(u) is smooth enough and hence {5n(u),v)jii is a linear functional 
w.r.t. V € U^{M.), which justifies identification of 5Yi{u) with an element of Z^#. Alternatively, the 
problem of finding the equilibrium configuration of atoms can formally be written as 

Mx^M, 



du{x) 
if we consider 11 as a function of finite number of variables u{x), x £ A^. 



For the equations (2.3) to have a unique solution, we must additionally require that the average 
of u is zero: 

{u)m = 0. (2.4) 



The equilibrium equations (2.3) together with the additional condition (2.4) can be written in 



variational form: find u G Wpor(A^) such that 

{5E{u),v)m = F{v) yvGUperiM) (2.5a) 

{u)m = 0, (2.5b) 

where the functional derivative 6E : Upcr{-M.) -^ Z^pcr(-A4) is computed as 

{6E{u),v)m= (y^V^r{D'R..u),Drv) , (2.6) 

and V^j.{Dfi^u) denotes, efi'ectively, the gradient of a scalar function V^ w.r.t. its vector-valued 
variable Dr-u (note the difi'erence with V^^fs introduced in ( |2.2[ )). Here and in what follows, with a 
slight abuse of notations, we keep the sign of summation over r E 7^^ inside the triangular brackets 
of the scalar product. 

2.2 Example: a Simplified Model 

The following simplified model will be useful in illustrating the numerical methods presented in 
this paper. 

Assume one space dimension, d = 1; the domain fi = (0, 1], the shift vectors in fixed configura- 
tion 

the multilattice 

m— 1 



M=[j{ez + e^)nn = ^znn, 



a=0 



Kl K2 km kl K-m-i K-m 

Figure 1: Illustration of a simplified atomistic model 



and the basic lattice C = e'Lr\Q,. We further assume TZ = {— } (nearest neighbor interaction only) 
and consider the "linear spring model" with the atomistic potential 

V,{DrU-x)=i^,{x)^-^^, (2.8) 

with r = —. Such system can be interpreted as a system of masses located at positions x + u and 
connected with ideal springs with spring constants k = 'ip)^{x)/e, as illustrated in Figure [I] 
The equilibrium equation then becomes 

{tp^DrU,DrV)M = {f,v)M (2-9) 

3 Quasicontinuum (QC) Method 

Traditionally, numerical methods such as the finite element method (FEM) are applied to continuum 
equations which can then be solved on a computer. The characteristic feature of atomistic models 
we are discussing in the paper is their discreteness, with a number of degrees of freedom often too 
large to keep track of each individual atom. Therefore, similarly to FEM, the ideas of reducing 
the number of degrees of freedom are used for atomistic models as well. The difference is now that 
the reduction is done from a large but finite number of degrees of freedom to a smaller number 
of degrees of freedom. The QC method is a representative of such methods. We first present its 
simple-lattice version. The QC method consists in reducing the number of degrees of freedom of 
the atomistic system by choosing a coarse mesh of nodal atoms and assuming that the positions of 
the other atoms can be reconstructed by a linear interpolation. 

It should be noted that we discuss here only the local version of QC which is equivalent to 
applying FEM to the Cauchy-Born continuum model of elasticity. We are not considering coupling 
the continuum and discrete models in this paper. 

3.1 Notations 

Assume a partition Th of the domain fi into simplicial elements T, which we will conveniently refer 
to as the mesh. Normally, #{Th) ^ #(>C) (recall that by #(•) we denote the number of elements 
in a set). By \T\ we denote the Lebesgue measure of T. The QC solution will be denoted as u . 
The space of piecewise linear discrete vector-functions is denoted as 

<er = {n' e {w^;,rmy ■■ At G Pi{T) VT g %}, (3.1) 

and the space of piecewise constant vector-functions as 

Qper = {?' e {L^,,mY : q'^lT G Po(T) VT G %} . 



3.2 QC for simple lattice 

In this (and only this) subsection we make the simple lattice assumption. That is, we assume that 
m = 1 and hence M. = C. In particular, in this subsection we write Ve{Dfiu;x) = V{Diiu) and 
TZeix) = TZ as they no longer depend on x. 

The QC method [39j aims at finding a minimizer of 

Uiu'') = {V{Dnu^))^-F{u^) 

in ZYper- Minimizing Il{u'^) in Up^^ indeed reduces the number of degrees of freedom of the sys- 
tem from 0{#{M)) to 0{#{Th)) (recah that #{%) < #{M)). However, one must stih spend 
0(#(A^)) operations to compute the effective forces on the reduced degrees of freedom. In order 
to have an efficient numerical method (i.e., a method with 0(#(7/i)) operations) one introduces an 
approximation to Ii{u^) which is called the local QC method [3H] (hereinafter referred to as the QC 
method). 

The local QC method first approximates DrU^ with Vr-u^ within each T (hence the name of the 
method: the nonlocal finite difference DrU is approximated with the "local" directional derivative 
VrU^)- Then for each x £ T one has 

V{Dnu^) « V{Vnu^) = W{Vu^\t), 

where M^(F) := V{FTZ) is the Cauchy-Born energy density associated with a displacement gradient 
F (see (A.4|) to obtain the precise definition of FT?.). Second, the local QC method changes the sum 



over X G £ effectively to integration over Q, i.e., 

E'i%u''):= f W{Vu'')dx= ^ |r|^(Vn'*|r). 

The variational formulation of the QC method is thus 



I ^5W{VrU^):Vrv'^dx = F'^{v^) MvdU^^,, (3.2) 



where 5W denotes the derivative of W and the semicolon denotes the inner product of matrices. 
For uniqueness of solution we should restrict ourselves to u G Uh. 

3.3 Multilattice QC 

Approximating the exact minimizer of 11 (u) with a piecewise linear u £ ^ner may be accurate 
enough for the case when the interatomic interaction V^(»,a;) varies smoothly with x (more pre- 
cisely, if the mesh Th resolves well the variations in V^(»,x)). However, for materials with multi- 
lattice structure (examples of such materials were given in the introduction) the piecewise linear 
approximation of the displacement u cannot be accurate. 

In this subsection we present the Multilattice QC (MQC) method first introduced in [30] which 
is designed to handle the multilattice microstructure. 

Define the space of QC deformations of the multilattice 

m—l 
Uh^<i = [yh+Yl <i>a ■■ U^ e <r, ll e Q^er, « = 1, • • • , m - l} , (3.3) 

a=l 



where g^ are the deformed shift vectors (recall that pa are the undeformed shift vectors) and 

Wa\c+epp = Sal3 («, /3 = 0, . . . , m - 1), (3.4) 

are the associated basis functions, with (5^^ denoting the Kronecker delta. For a more detailed 
introduction of the space of QC deformations, refer to [1]. In each element T G 7/i we have m — 1 
nonzero shift vectors q^ and we set Qq := 0. We denote 

Next, form the interaction energy E{u) with u = u^ + "^^Zi Qa'^a £ l^'^''^'- 

m—l 

E{u)= E[v^+Y,q^^Wo) 

a=l 

m—l 

K(^7^.w(n^x) + Y, Drql{x)w^{x))- 

a=l 
^ m—l m— i 

— XI ^4^n,{x+epp)[u^ix + epp) + 'Yq'^{x + €pf})wa{x + epf})j;x + epi3 

^ m—l m—l 

(—Yl ^^ADTi.,p^^{^ + ^Pt^) + Yl ^T^e,p(la{x + eP/3)Wa(eP/3) 

\ lit ^ •' ' XkizLf 

13=0 a=l 



where we used periodicity of V^ (see ( |2.2[ )) and Wa (which follows directly from the definitions of 
Wa and Ai). Similarly to the simple-lattice QC, we perform a local quasicontinuum approximation 
which consists in: (i) assuming that the energy associated with each T depends on the displacement 
gradient and shift vectors only in T, and (ii) changing the summation over x € C to the integration 
over fi: 

„ ^ 771-1 m—l 

E{u) - -Yl ye,JVn^,,u'' + Y <laDTi^,,Wa{epp))dx 

/o iif^ — ^ — ^ 

•^^^ /3=0 a=l 

^ m—l m—l 

= E 1^1 - E K,/3((Vn'^|T)7^.,^ + Y {q'a\T)D^^^^^w.{epp)) =: E'-^%u\c^'), 

T&Th /3=0 a=l 

where we used the identity Vt^^ a^^lr = (VM^|T)'7^e,/3, cf. (A. 4). 

Remark 3.1. The expression for E^'^'^ {u , {q^}) can be further simplified by denoting the species 
of atoms e(3 + TZ as A^^/s (formally A^^p '■= {ap^r)r£n^ ,3 where ap^r G {0, . . . , rn, — 1} is defined so 
that Pap r G p/3 + r + Zj. Then the sum in E^'^^ can be simplified as the difference between the shift 
vectors of interacting atoms: 

m—l , m—l \ 

Y {'ia\T)Dn,^Wa{epp) = I Y {Qa\T)DrWa{epi3)] 

a=l ^ a=l ' r&n^^ft 



(«,Jt)K(6p,,,J -u;.(6p,))) = [{q^ajr) - (^JIt)) 



r■e7^,,,3 \- p.r ' - ^ '/ren,,f) 



This yields 



m—l 



E-^^{uMq^j)=Y: I^I;^E^^./^((^^..^' + 4,,-'?')It)- 



(3.5) 



TGTh 



/3=0 



In the next step, the shift vectors qa are ehniinated from (3.3) by requiring the variation of 
^mqc^^/i^ l^^l^ w.r.t. Qy in each triangle be zero: 



m— 1 



m—l 



iEEK,,r{i'^AT)n,,+ Y.i^'\T)Dn. 



0=0 r&lp 



a=l 



\DrW^{epp) = 
(7 = l,2,...,m-l). 



(3.6) 



The equations (3.6) form a system of m — 1 equations for m — l unknowns {qa)^Zi in each T, from 
where under some stability assumptions (see, e.g., [20]) it is possible to determine the shift vectors 
qa depending (as a rule, nonlinearly) only on the displacement gradient: 

q^|r = q(Vn^|T). 

Note that the function q(F) does not depend on T, unless different periodic materials are considered 
in different elements T. 

We now form a QC energy with q^ eliminated: 



The QC equation of equilibrium now reads: find v!^ G ZYp^j. such that 



(3.7) 



{5E'^''^{u^),v^)n = F^{v^) ^v^ G U. 



h 
per' 



The function u gives a macroscopic displacement of the material, and one needs to compute 
u^ + '^^^Zi Qa'^a for the microstructure. We note that since qa were found from letting the variation 
of E^'^'^{u ,qa) w.r.t. qa be zero, we have 



5^mqc(^/.) = 5,J?'"^'(n\q(Vu'^)) 



(3.8) 



Remark 3.2. Instead of eliminating q = q(Vn ), one could also look for a critical point (or a 
minimizer) of the energy £'™^^(u^,q'*) w.r.t. both u^ and q^. In fS^ the latter approach is reported 
to be more efficient. 

3.4 Application of QC to the Simplified Model 



We illustrate an application of QC to the simplified ID model (2.8) for two species of atoms (i.e., 
m = 2), ^,(0) = ^1,^,(1) =i>2. 

If we approximate the exact solution with u G ZYp^j. as it is done in the simple-lattice QC (cf. 



Section 3.2) then we will find the approximate energy 



Em- 



T£Th 



V'l ^ ^ + 4^2 



Em 

Ten 



■01 + V'2 (V^n'^) 



h\2 



10 




k2 



Figure 2: Illustration of a 2D model problem with heterogeneous interaction. 



Here ip'^ = ^^^^^ is the wrong effective spring constant, since if the two springs in series are replaced 
with two identical springs with the effective spring constant ipo then V'o = .hfitii^ {see, e.g., [H]). 



V'l+V'2 



If instead we allow for nonzero shift vector gi then the corresponding MQC energy (3.5) is 

(V^n'* + g^)2 (Vru'^-q'^y 



E^^-{u\q'l)= Y, \T\\ 



T&Th 



Ipl 



+ ^2 



with ?" = 2- '^^^ strong form of (3.6) in this case can be obtained by differentiating the above 
expression w.r.t. ql in each T: 



V^i((V,n'' + q'i)\T) - V2((V,n'^ - gHIr) = 0, 



from where we find 






?llT 



Substituting this back to the MQC energy (cf. (3.7)) yields 



^mqc^^/, 



Em- 



T&Th 
Ten 



■01 



1 



2i^2 



2V01 + V2 



(V.^'^|t))' + 02J 



2tpi , ft, -^2 



2VV'i + V2 



2ViV'2 (V^ii'^i-^^ 



V'l + V'2 2 

where the effective spring constant is now computed correctly. 

4 Homogenization of Atomistic Media 

We now present another coarse graining strategy based on homogenization. We derive below the 
homogenized model of the atomistic material which will be the basis for formulating and analyzing a 
quasicontinuum method for multilattices. We will first present the "classical" approach of upscaling 



the atomistic equations to a continuum model and in Section 4.3 we will outline a strategy where 



the multilattice atomistic model is upscaled to the simple lattice model. 



11 



We note that the first works concerned with upscahng atomistic equations are [11^ [121 1231 • In 
the present section we derive the upscaled equations for a general model of interaction in many 
dimensions as opposed to the pairwise interaction in ID assumed in the upscaled equations |1H 
[T2| ^^ . The upscaled equations are derived using a formal asymptotic expansion. Rigorous error 
bounds for the homogenized equations can be found in the preprint [4J for the case of the ID nearest 
neighbor interaction and in [3] for the case of a ID finite-range nonlinear interaction. 

4.1 Asymptotic expansion 

In order to take into account the local variation of the atomistic interaction we think of the dis- 
placement as depending on a fast and a slow scale u{x) ~ u{x,x/e). We define x E M'', the macro 
("slow") variable, and y = x/e € Ij + V, the micro ("fast") variable, and consider a series of 
functions n" : M"^ x (Z'^ + V) -^W^ indexed by n = 0, 1, 2 . . . As we consider the local structure and 
interaction to be periodic, we assume that the functions u" are P-periodic in the fast variable, i.e., 
they satisfy for all {x,y) E Q x V 

u^{x,y + j) = vJ'{x,y), VjeZ'' 
while the behavior w.r.t. x is similar as considered in the previous sections 

n"(x + i,y) = n"(x,y), Vi G Z'^. 
We then consider the asymptotic expansion 

n(x)~ (n°(x) + en^(x,y) + eV(2;,y) + ...)|^^^/^ ^x e M. (4.1) 

Notice that we directly assume that the homogenized solution, u^, does not depend on y. 



We now proceed as in the "classical homogenization" [HI El ES] and plug the ansatz (4.1) into 



(2.5): 






y=x/el M 

where the test functions v = v{x,y) are continuous and smooth in x G and discrete in y £ V. 



Here we used the relation (A. 3) to expand the full derivative D^ through partial derivatives Dx^r 



Dy^r, and the translation operator Ty^r, and used the collection-of-derivatives notation D-ji (see 



Appendix A. 3 for more details). 



We then extend the equation on the entire Ai x V: 

( Y, Vl.iDx^n.u'^ + eDx^nTy,n.u^ + Dy^n^ + • • • ) 
veil, (4.2) 

, Dx,rTy,rV + e^'^Dy-rv) = (/, v)mxV- 

I MxV 

We now expand this equation in powers of e. For that, we use the approximation Dx,r ~ ^x,r 
(i.e., we essentially use Taylor series to expand D^^r), and the notations Ve{»;x) = V{»;y) and 



12 



TZe{x) = Tl{y), and change a sum over Ai to an integral: 






nx-p 

where {•)n,v is a short-hand for J^(»)-pdx. 



We first cohect the 0(e ^) terms in (4.3): 



V V; (V,,7eu° + Dy^Tzu^) , Dy,.v) = 0. 
reTe 

As usual in homogenization we write the solution of this equation (of course, equipped with the zero- 
average boundary conditions) as u^{x, y) = x{^ xU^ {x); y) + u^{x), where x = x(F; y) '■ W^^'^ xV ^ 
W^ solves 



findx(F,.)e^/#(7')s.t. (^V;{Fn + Dy,nx{n),Dy,r<r) =0 ya€U#{V). (4.4) 






To obtain the equation for the homogenized solution u^{x), we collect the 0{e^) terms in (4.3) 
and use the test function i; of x only: 

( Yl K'(Vx,7^^i° + Dy^Tzu^),V^,rv')^^ = {f,v)nxv- 

This leads to the homogenized equation 

(<^$0(V,nO),V,.^)n = (/,^)n, (4.5) 

or equivalently in the strong form —Vx ■ S^^{Vxu'^) = f{x), where (5$° : M'^^'' —?- M'^^'^ satisfies 

6^\n = ( Yl K'(F7^ + Dy,nx{n) r") ^^- (4.6) 

reTl(y) ^^ 

Thus, we obtained the equation for the homogenized displacement u^ with the homogenized tensor 
6^ . For the equation (4.5) to have a unique solution, it has to be supplemented with boundary 
conditions, for instance by requiring that u^ is periodic and has zero average. 

As an illustrative example, in the case of ID with a pair interaction potential we can write 
V{DTi(^y-^u\y) = Ylrell(v) ^r{DrU;y) (cf. the Lennard-Jones potential in (2.1)), consequently. 



reTl(y) ^^ 

where the derivative of $r is with respect to F. 

Remark 4.1. In the above formal arguments we assumed, for simplicity, that the external force 
f is a continuous function of x and moreover does not depend on y, i.e., f = f{x), x € Q. We 
emphasize that oscillatory external forces could also be considered. The homogenized equation would 
then depend on a proper average of the external forces. The assumption that f is a function of the 
continuous variable x can later be relaxed once the homogenized equations are discretized on a finite 
element mesh. 
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4.1.1 Underlying Homogenized Energy 

It is useful to highlight one more feature of the homogenized equations. 

Proposition 4.1. The function 6^ (f) defined by (4.6) is the derivative of the following function 



$0(F):=(y(F7^ + Z),,7^x(F)))^,^ 
w.r.t. F. 
Proof. Compute the variation of ( |4.7[ ) w.r.t. F: 

6<^\f) : G = ( j; y; (F7^ + Dy^nxin) ' (Gr + Dy^nSxi^ ■ G) 
re7^ 



yer 



(4.7) 



(4.8) 



Since dx{^) :G G U^{V), the second term in (4.8) drops due to (|4.4[) and we have 



5cl>0(F):G = ( Y, V;ifn + Dy,nx{^)) ■ Gr) ^ = ( J^ y;(F7^ + A,,7^x(F))rT :G 



re7^ 



/j/eP 



r•e7^ 



j/eP 



which is consistent with (4.6). 



D 



Corollary 4.2. The equations (4.5) can be written as 

{6E\u^),v)n = {f,v)n, 

where 

^O(nO) := [ $0(Vn°)dx. 
Jn 



(4.9) 



The fact that the homogenized equations have an underlying energy may be important in some 
applications where, for instance, one chooses to use nonlinear conjugate gradient algorithms or 
needs to check for stability of numerical solutions. 

4.2 Application of Homogenization to the Simplified Model 



For the simplified model described in Section 2.2, the cell problem (4.4) in the strong form reads 



Dy^.ri^ (Fr + Dy^rXiF))) = 0. 
where r = — and ip{y) := ipei^y)- Its solution can be written as 



Dy,rXiF) 



C 



Fr. 



with C = Fr {l/ip)-p, and V given by (2.7). The homogenized energy density is therefore 

$0 = (^^^Fr + Dy,rX)\ = {^Fr +^-Frf)^^^-- ^' 
which yields the homogenized energy 



(iM 



V 



1 (Fr)^ 



(4.10) 



(4.11) 



i?0(n°) 



(iMp' 



(Vrti' 



0^2 



dx, 



with the correct form of the effective spring constant ij)^ = {\/'ip).p . 

We emphasize that this procedure and the obtained results are well-known for PDEs [8", Chap. 



1], and are in agreement with the equations obtained using the multilattice CB rule (Section 3.3). 
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4.3 Discrete Homogenization 

Homogenization effectively upscales the discrete problem with microstructure to a continuous prob- 
lem of nonlinear elasticity. One can, however consider another approach where the upscaled model 
is a discrete model with no microstructure. 



For that, one assumes that x is discrete, x G £, and approximates D^^r* ~ {Dx*)r in (4.2) 



where Dxu{x) G M is the discrete gradient of n G Z//per(>C) at the point x £ C defined as 
{Dxu{x))ek = Dx^e^u^x), k = 1, . . . ,d. Following the above procedure of asymptotic expansion one 
gets the upscaled equation 

{S^°{Dxu''),Dxv)c = {f,v)c yv£U#{C), (4.12) 



where 5$ is defined by ( |4.6[ ), the same equation as for the continuum homogenization. 

Thus, we obtained the equation for the homogenized displacement u^ with the homogenized 
tensor 6^ . The homogenized tensor 5$ no longer depends on the fast variable y and therefore we 



can apply the standard QC method to the homogenized equation (4.12). Again, to ensure a unique 
solution for (4.12) we search for u^ G IA^{C) (see the preprint |4j for more details). 



One of the advantages of the discrete homogenization is that it is not required to assume a 
continuous force /. Another advantage is that the discrete homogenization can potentially be 
helpful in deriving effective interaction potentials, for instance by averaging over the temperature- 
related oscillations of atoms around their equilibrium positions. 

5 Homogenized QC 

We formulate a numerical macro-to-micro method for treating multilattices, which we call the 
homogenized quasicontinuum method (HQC). We introduce HQC in the framework of numerical 
homogenization. For the case of materials with known periodic structure (i.e., crystalline materials) 
the HQC method will be shown to be equivalent to applying finite elements to the homogenized 



equations (see Theorem 6.1). 



Nevertheless, we argue that HQC can be applied to non-crystalline materials and to time- 
dependent zero-temperature and, possibly, finite-temperature problems. Indeed, in Section [7] we 
give an application of HQC to a stochastic material and in Section [8] we present an application of 
HQC to a ID time-dependent zero-temperature evolution. In addition, the HQC serves a convenient 
framework for the error analysis |3l U] . 

We present the HQC algorithm assuming that the microstructure is a function of the macro- 
scopic displacement. A reformulation analogous to the concurrent coupling of [37J is also possible 



(cf. also Remark 3.2) 



5.1 HQC Method 

The method will be presented using macro-to-micro framework as used in some numerical homoge- 
nization procedures [U [IBl [261 E2l SI] • We present the method for the case when the external force 
f = fe may be microstructure-dependent. 
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5.1.1 Macroscopic afRne displacement 

We again assume a partition Th of the domain il into simplicial elements T, recall the definition of 
the space U^^^., (3.1), and introduce its subspace of zero-mean functions U41 C Z^per- 



5.1.2 Sampling Domains 

We choose a representative position x^'^ G £ and a sampling domain S^^ := x^'' + eV associated 
with each T £ Th- The sampling domain is normally chosen inside T (the mesh can be highly refined 
in certain regions and therefore some sampling domains S^^ may be bigger than T). However, in the 
problems of coupling fully atomistic and coarse-grained equations, the mesh may be non-uniform 
and some T £ Th may be too small to accommodate a sampling domain fully inside T. 

The sampling domains have the associated operator of averaging over the sampling domain. 



(•)^g5.rcp and the functional space U^i^S^'^) = U^{eV) (see (A.l) for the precise definition). 



5.1.3 Energy and Macro Nonlinear Form 

Define the atomistic interaction energy of the HQC method 



E^^^iu') := ^ \T\{V,iDn,RTiu')))^^^..,, (5.1) 



-'T 

TeTh 



Qrep 

iJrp 



where Rt{u ), defined by (5.3), is the microfunction constrained by u in the sampling domain 



The functional derivative of the above energy reads 



{SE^^%u'^),v% = Y^ \T\( Y^ V^^,{Dn^RT{u^)),n,5RT{u^)v^) . (5.2) 

where 5Rt{u^) is the functional derivative of the reconstruction Rt{u^) defined below. 



5.1.4 Microproblem 

Given a function u^ G ^peri Rt{u^) is a function such that Rt{u^) — Uy^ G h(^{S!^^) and 

YvlriDn.RHu'')), Drs) =0 V5GZY#(5^^P), (5.3) 



where u[^j^ is an affine extrapolation of u^\t over the entire M"'. If S!^^ C T then u^^^ can be 
substituted with u^. 

Remark 5.1. When modeling essentially nonlinear phenomena (e.g., martensite-austenite phase 
transformation), one should require that the micro structure corresponds to a stable equilibrium. 
That is, one should require, in addition to (5.3), that w = Rt{u^) — u{jn G U^{S^^) is a local 



minimum of {VeiDnSuL + w))) ^^s'^^p g^ p. 238]. 

Remark 5.2. In the case of linear interaction, the reconstruction Rt is a linear function and hence 
5Rt{u^)v^ = Rt{v^), which makes the derivative of the HQC energy (|5.2[) take the form 






'^ ' 
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Remark 5.3. The functional derivative of the HQC energy (5.2) can equivalently he written as 



{SE'^^'^iu'),v% = E l^l( E y:ADn.RTiu')),iVrv'\T)) _, (5.4) 



TeTh ran, ^ ^ 



by noting that Dr5RT{u^)v^ = ^r^^un + [Dr5RT{u^)v^) - DrV^^^^, that 



m 



view of (5.3), and that Dr-Vy^^ = Vr-w on each T. Here we used the fact that 5Rt{u )v — fjj^ G 



repN 



U^{S^^) which follows from taking the functional derivative of Rt{u^) — Uy^ € U^{S!^ 

5.1.5 Reconstruction 

The functions Rt{u ) describe the microstructure of the solution inside each 5'^''. One can re- 
construct the solution describing the microstructure, u^''^, from the homogenized solution u^ by 
combining Rt{u ) into a single function defined on entire atomistic lattice M: 

u^^^{x) = Rt{u^){x) (xGTnM). (5.5) 

That is, we effectively extend Rt{u^) periodically on each T. It should be noted that (5.5) does 



not uniquely determine u '^(x) if x G dT for some T ^Th- 

5.1.6 Variational Problem 

We define the homogenized quasicontinuum approximation as the solution u^ E Uh of 

{5E^^^{u'^),v'^)n = F^'i^v'') ^v^ £ Ul (5.6) 

where 

^hqc(^h)^ Y^\T\{f,,v^),^s--- (5.7) 

TeTh 
If the external force is smooth, it could instead be evaluated for a single representative atom. 

In the case of linear nearest-neighbor ID interaction it can be shown that (5.7) is well-posed 
and that the homogenized quasicontinuum solution u^ approximates the solution u of the original 
equations only in the L2-norm. To get a good approximation in the //^-norm, the reconstructed 
solution n^'^ should instead be considered. We will report the analysis for the nonlinear case in a 
separate paper (see the preprint [U Theorems 4 and 5] for the analysis of a linear model) . 

5.2 HQC Algorithm 



The problem (5.6) is nonlinear and its practical implementation is usually done by the Newton's 



method. We briefly sketch below an algorithm for solving (5.6) 



For the Newton's method we need to compute the second derivative of the energy (5.1): 



{5^E'''\u'')w\v^)n= E |r|/ Y. ye[r,p{Dn.RT{u^))Dp5RT{u'^)w\ Dr6RT{n'')v^ 



T 



(5.8) 
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5.2.1 Newton's Iterations for the Macroproblem 

The algorithm based on the Newton's method consists in choosing an initial guess u '^^> G [/i and 
performing iterations 

(52i^'^^'(tx'^'("))(n^'("+i)-n'^'(")),t;'^>^ = Fi^<^^(z;'') ^v^^eU^, (5.9) 

until u^'("'~^^) becomes close to vf^'^"^' in a chosen norm. 

To solve the linear system (5.9) for vr^\^~^^i — u A'^) g ]j'l^ we choose a nodal basis w^ (1 < ^ < 



K) of Wpgj,. One way to satisfy the condition {u )q, = would be to perform all the computations 
with one basis function eliminated (e.g., to consider w^ for 2 < k < K), and post-process the final 
solution as u^ — {u^)n. 



The stifi'ness matrix of the system (5.9) will thus be 

and the load vector will be 

b^ = F'^^^(u;i). 

As given by the formula (5.8) we need to compute the solution of microproblem Rt(u^'^"') on each 



sampling domain S^^ as well as its derivative 5Rt{u^'^"'^)w^. 
5.2.2 Solution of the Microproblem 



The microproblem (5.3) can also be solved with Newton's method. For that, in each T one needs 



to choose an initial guess n*^"'*^' to Rt{u^'^^'), for instance u^^'^'{x) := n'*'*-"'(x) and solve 

v5GZY#(5n, 

with respect to n^"'^"^^-* constrained by u("''^+^) — Uj^^ G Z^#(-S'^''), until the difi^erence between 
y^{n,u+i) g^j^i^ y(n,u) ^g gj-Q^H vn a chosen norm. 

After that, we can compute 5RtWi = 6Rt{u '^^')wi by solving 

Y, y::rADn.u^"'''^)D,i6RTwh, Drs) =0 VsGW#(5n (5.10) 

constrained by 6Rtw^ — (iff )iin S U:ff,{S^^). Notice that the gradient of all but d+ 1 basis functions 
Dr{'Wi)iin inside T are zero, which implies that we essentially need to solve the problem (5.10) d + 1 
times. 

Also observe that when computing 6Rt{u '^^')wi, we need to invert the same linear operator 
as in the final Newton's iteration, which allows for some additional optimization. 
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5.2.3 Possible Modifications of the Algorithm 

First, notice that when solving for u '^^~^^> we could linearize the problem on the previous iteration 
y^h,{n) ^ In that case we would have linear cell problems and thus we would need only outer Newton's 
iteration, but it would be required to keep the values of the micro-solution Rt{u '^^') from the 
previous iteration. We notice however that for a practical implementation of the above algorithm 
it may also be required to keep the values of the micro-solution: one needs these values to initialize 
the inner Newton iterations; depending on the initial guess for the microproblem the iterations may 
converge to a wrong microstructure. 



Another modification could be to compute the contribution of the external force /^ in (5.7) for 
a single atom in the case of no oscillations in Z^. 

In the case of linear interaction, the algorithm becomes simpler: one does not need to do Newton 



iterations. Nevertheless, the algorithm in Section 5.2 is applicable to the linear problem where it 
converges in just one iteration. 

6 Equivalence of Numerical Methods for Multilattices 

In this section we show the equivalence of three different methods for computing equilibrium of 
multilattice crystals, namely (1) the proposed HQC method, (2) finite element discretization of 
continuum homogenization, and (3) multilattice QC. We assume that in all four methods the 
external force / is approximated by the same function f^. 

Below we specify the four methods that we compare. It should be noted that given the macro- 
scopic displacement u we cannot guarantee uniqueness of the energy as there may be several 
solutions to the micro-problems corresponding to different phases of a multilattice crystal. To ad- 
dress such non-uniqueness, we treat the energy and all the respective microfunctions as multi- valued 
(i.e., set-valued) functions of u^ G Uht and prove that these multi- valued functions coincide on each 



Method 1. (HQC) We define the energy of the HQC method, £;^^^(u'*), by (|5l|, where Rriv!') 



is the set of all solutions of (5.3). 



Method 2. (FEM for homogenized equations) The energy of FEM discretization of the ho- 



mogenized energy is E^{u^) given by (4.9), and (|4.7[), and x = xi^'^v) is a set of all solutions 



of (4.4). 



Method 3. (Multilattice QC) E'^'^''{u^) is defined by (|3^ and q'^ = q(VM'') by K^. 



Theorem 6.1. E'^'i''{v!') = E^{u^) = E'^'i'^{u^) for any u^ e UJ^;^,. 

Proof. Part 1, E '^'^{u ) = E^{u ). First, we show that the micro-functions of Methods 1 and 2, 
Rt and x, are related through 

{Rt{u'')){x) = <(x) + 6x(Vn^|T; f ). (6.1) 

Indeed, denote F = Vu^\t and compute DrRriu^)'- 

DrRrin'') = D.u^ + eZ?,x(F; f ) = Fr + Z)j,,.x(F; f ). (6.2) 



19 



The following calculation shows that the left-hand sides of (5.3) and (4.4) coincide up to a factor 



Y,ye,r{Dn.{x)RT{u'')-,x),DrS{x) 

ren 



x&S^ 



Y, V;{DT^^y)RTiu^); y), e-^Dy^rs{ey) 



rell 



yer 



«"'( Yl ^^(F7^ + Dy,n X(F; yy,y),Dy,ra{y) 



r•e7^ 



yev 



where we do the change of the independent variable y = -, and of the test function o"(y) = s(ey). 
Hence (6.1) indeed relates the set of solutions of (5.3) and (4.4) with F = Vu \t- 



The following straightforward calculation concludes the proof of E^^'^^u^) = E'^^ 

Ten 



Y, |T|(K((Vn'^|T)7^e + I),,7^.x(Vu^|T;f))L 



TeTh 



es^^p 



^ \T\{V{{Vu^\T)n + Dy^^^xi'^u^\T;y)))y^^, 

T&Th 

Y^ \T\^'^{Vu^\t) = E\u''), 
TeTh 



where we used (6.2) in the first step of this calculation. 

Part 2, E^'^'Hvr) = E^'^^{u^). The main component of the proof consists in fixing T £ Th and 
showing that q'^lx and Rt{u^) are related through 

(7^|t = f/(epa) - ^(0), a = 0,...,m-l, (6.3) 

where U := Rt{u'^) - u(j^ G U#{St^)- 

First, we define g^jr := U{epa) — U{0). Notice that due to eT^-periodicity of U, we can write 

m— 1 



U{x) = Y U{epa)Wa{x), 

subtracting the constant f7(0) and applying Dr yields 

m— 1 
DrU{x) = A.( - U{0) + Y, U{epa)Wa{x) 

0=0 

m—l 

= DrlY ^i^Pa) Wa{x 
a=l 
m—l 

= Dr Y(^a\T)Wa{x), 

Q = l 

where we used the identity ^^Zo ^a(^) = 1 for all a; G A^. 
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Figure 3: Equivalence of different methods 



We then substitute q^\T = q^\T into (3.6). The argument of Fe in (3.6) can be written as 



m—l 



m—1 



{Vu^\T)n,,^ + Yl {q^\T)Dn,^^Wa{epp) = Dn^^Ju^ + Y^^qak 



a=l 



Wr 



a=l 



(6.4) 



Dn,,, (< + U{ep^)) = Dt^^^,)Rt{u^){ 



\x=ep0 



and therefore, upon noticing that summations over x = epfs and over x G S"^^ coincide for the 
eP-periodic functions, we conclude that the left-hand sides of (3.6) and (5.3) coincide when s{x) is 
chosen as s{x) = w^{x) — {w^{x))xee'Pi j = I ■ ■ ■ ,m — 1 (then Dr-s = Dr-Wj). This proves that ^^|t 
satisfies (3.6), i.e., g^|r ^ QoIt = U{epa) — U{0) (in the sense of the sets of solutions). 

To show the converse we notice that (5.3) holds with the function s{x) = w^{x), j = 1 . . . ,m — l 
and, obviously, with the function s(x) = 1. These functions form a basis oiUperi'P) = l^pcriS^"^), 
therefore (5.3) holds with any s G h(^{S^^) C Uper{S!^'^). Hence, g^|r C U{epa) — U{0), which 



concludes the proof of (6.3) 



The stated identity E'^'i''{u'') = E'^'i''{u^) follows directly from K^ 



D 



Remark 6.1. One can consider yet another approach to coarse- graining multilattices, namely 



consider the discretely homogenized method (4.12) and apply the standard QC method (see Section 
3.2) to it. As a result we will obtain energy of {^^{'Vu^))q which obviously coincides with the 



energy of FEM applied to the continuously homogenized equations. 



As a corollary of Theorem |6.1| and Remark 6.1, the solutions corresponding to the different 



methods considered, being critical points of the energy, also coincide (of course, provided that the 



external force is treated in the same way for these methods). The Theorem 6.1 is graphically 
summarized in Fig. |3J 
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Figure 4: Stochastic atomistic model: An illustration of the model (left) and an exact solution for 
32 X 32 atoms (right). 



7 Application of HQC to Stochastic Materials 

The HQC method can readily be generalized for non-crystalline materials such as glasses or complex 
metallic alloys. For that, lacking the period of the microstructure V, one only needs to take S'^'' 
large enough to accurately represent the material's microstructure. In this section we present an 
example of such computation. 

Additionally to taking 5^^ large enough, one could also average over an ensemble of samples 
of different microstructures for a given macroscopic displacement gradient Vu \t in each element 
T; however, we do not pursue this in the present work. We refer to |101 [27] and references therein 
for theoretical studies of stochastic homogenization of lattice energies. 

We take an atomistic system of 2048 x 2048 atoms. The atomistic bonds are chosen to have 
quadratic interaction energy. 



^(^) = (5^^V'e,r(x)(Z?.n) 



reTe 



x&M 



with TZ = {(1, 0), (0, 1), (1, 1), (— 1, 1)}, as illustrated on Fig. 4(a) The bonds' strengths Ve,r are 
randomly generated with a uniform distribution between 0.5 and 10 for r = (1,0) and r = (0,1) 
(i.e., vertical and horizontal bonds) and between 0.1 and 5 for r = (1,1) and r = (—1,1) (i.e., 
diagonal bonds). The external force is chosen as 



/(xi,X2) = i0e-^°'^("^i)'-™'^("^2)2 



sin(27rxi) 
sin(27rx2) 



/, 



where / is determined so that the average of / is zero. The equilibrium configuration for a system 
with 32 X 32 atoms is illustrated on Fig. |4(b)[ 

We then apply the HQC algorithm to that system. For the microproblem we choose a subsystem 
of A'rep X -/Vrep atoms and precompute the effective elasticity tensor. We then compute the HQC 
solution and compare it to the exact solution of the problem. A structured triangular uniform mesh 
with right-angled triangular elements with the leg size h = j, o, . . . is used. 
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Figure 5: Dependence of relative error of computing the energy with HQC and a straightforward 
appHcation of the Cauchy-Born rule. The squares and diamonds correspond to N^ 
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2045 



I.e., 



when the microproblem coincides with the entire system). A second-order convergence of HQC is 
observed whereas the Cauchy-Born rule solution does not converge. The dot-dashed and dotted 
curves are \E!^'^'^ — E\ for A^rcp = 128 and A'rep = 16 respectively. 

For comparison, we also produce the results of calculation with (a straightforward application 
of) the Cauchy-Born rule for computing the effective elasticity tensor; i.e., when atoms are not 
allowed to relax to equilibrium when an external displacement gradient F is applied. 

The relative errors of the interaction energy of HQC and CB solutions {Er^'^^ and S , respec- 
tively) as compared to the energy of the exact solution E^ are plotted in Fig. plfor different mesh 
size h and different sampling domain size A^rcp- A second-order convergence of HQC and absence of 
convergence of the solution computed according to the Cauchy-Born rule can be observed. One can 
also see that with iVrep = 128 (and even with A'rep = 16) one can get a rather accurate numerical 
solution. 

8 Application of HQC to Time-dependent Problems 

We apply the proposed HQC method to the ID zero-temperature evolution of a multilattice, de- 
scribed by the following equations 



u|j=o = u^ 
u\t=Q = 0. 



?.la) 
Lib) 
B.lc) 



where u = u(t,x) £ C^([0, T];Z^pcr(A^)) is the time-dependent displacement of atom x, vP 
vP{x) G Uper{-M.) is the initial displacement, M'^(x) = M(|) is the mass of atom x, ii 



di^' 



■u = gWii. We assume no external forces. 



One can, assuming no fast oscillations in time of the microstructure, perform the two-scale 
expansion procedure for the time-depend case (which closely follows the continuum case [8]) 

{M'^U, v)m = {5E\u),v)m Vw G ^/per(-M), 
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where E'^{u) is given by (4.9), and M" = {M)p and likewise formulate the macro-to-micro dis- 
cretization [21 [22] 

For the numerical test we take the same lattices as for the simplified model with m = 2 (see 



Section 2.2). The atoms interact with the Lennard- Jones potential (2.1) with 



is half-integer I 1.01 | is half-integer 

is integer, ' 1 0.99 f is integer, 



and the cut-off distance R = 3. The masses of atoms are 

, ^^ , , ( 1 - is half-integer 
M'{x) = < ' 

I 2 I is integer. 

The atomistic system contains j^{A4) = 2^^ atoms. 

The initial displacement is chosen in the following way: First, we compute an equilibrium 
displacement u, i.e., such that {5E{u),v)m = Vv G Uperi-M.)- Second, we compute an eigenvector 
of 5^E(u), ui, corresponding to the mode oscillating most slowly. Then, the initial displacement is 
taken to be n" = u + 0.01 u£)^]\ ^ ■ With such an initial displacement, the solution remains smooth 
(i.e., most of energy of the solution is contained in long wavelength modes) for times comparable 
to the oscillation period, and one can compare a QC approximation of the solution with the exact 
solution. 

We compare the reconstructed solution obtained by the HQC discretization in space with the 
reference solution obtained in the full atomistic computation. The reconstruction of the HQC 



solution is performed similarly as described in Section 5.1.5 The HQC discretization is performed 



on a sequence of meshes with h = | ,„,... . For the time integration, we use the Verlet method with 
the timestep r = ^h for the HQC solution and r = ^^je for the reference atomistic solution. We 
run the computation until T = ^, which corresponds to about a quarter of a period of oscillation 
of the solution. 

The errors in (the discrete analogs of) L°^([0, T]; L^(ri)) and L^{[0,T]; H^{Q)) are presented in 
Fig.pl One can clearly observe for relatively large h a second order convergence in the L^{Q)-noTioa 
and a first order convergence in the H^ {Q)-noTiR, and the convergence seems to stagnate as h is 
further reduced. 

9 Summary and Concluding Remarks 

We have considered the problem of equilibrium of multilattice crystalline materials and discussed 
the application of the (local) QC method [40j for such materials. We then have proposed a homog- 
enization framework and based on it proposed a numerical macro-to-micro method which we called 
HQC. We have shown that the four methods, namely the HQC method, the QC method applied to 
the discretely homogenized equations, the multilattice QC, and the finite element method applied 
to continuously homogenized equations, are equivalent. 

Despite equivalence of the methods for statics of multilattice, we argue that the homogeniza- 
tion framework developed in this paper has several advantages. First, it contributes to a better 
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Figure 6: Error of time-dependent solution in (the discrete analogs of) the norms L°°([0, T]; L^(il)) 
(left) and L'^{[0,T]; H^Q)) (right). 

understanding of the multilattice QC method and provides a link to the existing theory of homog- 
enization of PDEs. In particular, we have generalized and applied the HQC method to the case 
of random materials and to the unsteady case, numerically demonstrating convergence of the pro- 
posed numerical method. Second, the developed homogenization framework allows for application 
of analytical techniques available in the homogenization theory and thus seems most promising for 
convergence analysis of numerical methods for multilattices. We refer to our preprint [4J and an 
ongoing work [3] for an example of such analysis. We also note that the extension of the homog- 
enization technique proposed in this paper to atomistic materials at finite temperature is of high 
interest. 

A Notations 

In this appendix we gather the frequently used notations. 

A.l Function spaces 

For any finite set S C M , we define the discrete averaging (integration) operator {•)s by 

and sometimes, more verbosely, as {u(x))xes- Here i^{S) is the number of elements in the set S. 

We consider discrete periodic functions (e.g., displacements or external forces) with the periodic 
cell rj = (0, l]'' (d € N), and the lattice (being, actually, the discrete periodic cell) S C 0, {S = C,A4) 
containing a finite number of points: #(5*) < oo. The periodic extension of the lattice is denoted 
as Sper = S + Z'^. Such space of periodic functions is denoted as 

KpcviS) = {u : 5pcr ^ M : u{x + a) = u{x) Vx G 5, Va G Z'^}, (A.l) 
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and the space of periodic functions with zero average: 

U#{S) = {u^U^,,{S): {u)s = Q]. 

We do not have separate notations for scalar and vector-valued functions and explicitly state 
whether the function is scalar or vector-valued when it may cause ambiguity. 

Similarly to the discrete averaging, we also use continuum averaging notation (u)n := j^ u{x)dx, 
and for functions of two variables we write {v)sixS2 '■~ ((^)s'2)'? i where each Si {i = 1,2) can be 
either continuous or discrete. 

For vector- valued u = u{x) and v = v{x) we denote the pointwise scalar product as u ■ v (i.e., 
(u • v){x) = u{x) • v{x)) and the semi-inner product in Wpcr(/3) as 

{u, v)c = {u ■ v)c = -jrrw. XI ^(^) ■ ^(^)- 

(It is a proper inner product only in U#{C).) We similarly define the pointwise scalar product and 
the (semi-)inner product for functions of continuum variable and for functions of several continuum 
or discrete variables. 

A. 2 Operators 

For u : S ^M.'^ {S = C,A4) we introduce the finite difference Dx^rU 

Dx ru{x) := ^(^ + gO-^(^) (fo^ ^(^s, r G M'^ such that x + er € S). 

e 

In addition to differentiation operators, we define for u £ Wpcr(/^i), the translation operator TxU G 

^pcr(-Cl) 

Tx^ru{x) := u{x + er) (for x (^ S, r (^M. such that x + er G S*). 

The definitions of the discrete derivative and translation generalize to functions of two variables 
by considering the partial discrete derivative and translation operators, i.e., Dx^r^Tx^r applied to 
u{»,y) and Dy^r,Ty^r applied to u{x,»). 

In homogenization we consider "traces on diagonal" of function of two variables, v = v(x, -). 
For such functions we introduce full translation and full derivative operators Tr := Tx^rTy^r, D^ := 
\{Tr — I) SO that 



{TrU)\y^x=Tx;r\u\y^xy and {DrU)\y^x = Dx,r\u\y^xy (A. 2) 

The following relates the partial and the full translations and derivatives: 

Notice that the variables x and y are not symmetric in the definition of full derivative. If a 
function does not depend on y then the full derivative coincides with the derivative in x (likewise for 
the translation). Hence, for functions of x only, we sometimes omit the subscript x in the operators 

For continuous functions we denote Vu a gradient of u and VrU = (Vn) • r a directional 
derivative. For a vector-valued function u, the directional derivative, VrU is defined componentwise 
and the gradient S/u is a matrix such that S/rU = (Vu)r. 
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A. 3 Functions of Vector-indexed Variables 

We consider a general form of interaction, where an energy of each atom depends arbitrarily on 
relative displacements of all the nearby atoms. Namely, for the "interaction neighborhood" TZ = 
{ri, . . . , Tfc} we consider functions 

Since the interaction neighborhood may be different for different atoms (recall that we consider 
multilattices) and contain different number of neighbors fc, we index derivatives directly with r £ TZ. 
That is, we use the following notation for tuples a indexed with r £ TZ: 

(ar)re7e := ("ri, • • • , OrJ for TZ = {ri,...,rk} 

and define 

Dtzu := {Dru)reTZ, VTeu := (Vr^i)re7^• 

Thus, for the functions of 7^-indexed tuples we write 

V{Dtzu) := V {Dr^u, Dr^u, . . . , Dr^u) . 

The common algebraic operations on 7^-indexed tuples are taken componentwise, e.g.: 

Dtiu + Dtiv = (DrU + Drv)r(zn, F7^ = (Fr)re7e etc., (A. 4) 

which is fully analogous to the algebraic operations on /c-dimensional vectors. 
A partial derivative of V{D'jiu) w.r.t. DrU (r G TZ) is denoted as Vf^^D-jiu). 
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